#ifndef _SPEC_H_
#define _SPEC_H_

#include <fftw3.h>
#include <vector>
#include <complex>
#include <algorithm>
#include <functional>
#include <numeric>
#include <iostream>

#include "criterion.h"

template
<
    typename Compression
>
class OutBandEmission
{
public:

    const char* getName()
    {
        const static char name[] = "OutBandEmission"; 
        return name;
    }
    
    template
    <
        typename SignalSystem
    >
    double operator()(SignalSystem& system, double compression = Compression().value)
    {
        using namespace std;
        
        if(system.empty())
            return 0;
    
        vector<double> density(system[0].size());
    
        for(int i = 0; i < system.size(); ++i)
        {
            vector<complex<double>> signal(system[i].size(), 0);
            copy(system[i].begin(), system[i].end(), signal.begin());
            
            std::vector<std::complex<double>> spec = OutBandEmission::fft(signal);
        
            transform(
                spec.begin(), spec.end(),
                density.begin(), density.begin(),
                OutBandEmission::add
            );
        }
    
        transform(
            density.begin(), density.end(),
            density.begin(),
            bind2nd(divides<double>(), system.size())
        );
    
        double power =
            accumulate(
                density.begin(),
                density.end(),
                0.0
            );
        
        double outbandPower =
            accumulate(
                density.begin() + min((double)density.size(), density.size() * compression / 2.0),
                density.begin() + min((double)density.size(), density.size() * (1.0 - compression / 2.0)),
                0.0);
    
        return outbandPower / power;
    }

private:

    static std::vector<std::complex<double> > fft(const std::vector<std::complex<double> >& signal)
    {
        using namespace std;
        
        std::vector<complex<double> > result(signal.size());

        // Create plans
        fftw_plan plan = fftw_plan_dft_1d(
            min(signal.size(), result.size()),
            (fftw_complex*)(&signal[0]),
            (fftw_complex*)(&result[0]),
            FFTW_FORWARD,
            FFTW_ESTIMATE);

        // Forward DFT
        fftw_execute(plan);
    
        // Free memory
        fftw_destroy_plan(plan);
    
        return result;
    }

    static double add(std::complex<double> x, double y)
    {
        auto n = abs(x);
        return n * n + y;
    }
};

#endif // _SPEC_H_
